Repositorio de código

Información general sobre los datos

Durante el proceso de análisis, hemos encontrado que la base de datos correcta para realizar la anotación es la Ensembl en su versión 75 (2014), las versiones más recientes no contienen algunos de los genes que aparecen en el dataset. Por tanto utilizaremos esa para anotar los genes y calcular los GRanges para su ploteo posterior.

Selección aleatoria de los datos

En primer lugar vamos a seleccionar diez muestras aleatorias por cada grupo, primero seleccionaremos sobre el dataset de targets y luego cogeremos los conteos correspondientes a esas muestras.

# Only select and save the samples once, and set the sample seed just in case we have to resamples the same population
set.seed(28052020)
selected_data_file <- 'selected_data.Rdata'

if (!file.exists(selected_data_file)) {
  targets <- read_csv('../data/targets.csv') %>%
    group_by(Group) %>%
    sample_n(10) %>%
    arrange(Group) %>%
    mutate(rn = Sample_Name) %>%
    column_to_rownames('rn')

  counts <- read_delim('../data/counts.csv', delim = ';') %>%
    dplyr::select(X1, rownames(targets)) %>%
    column_to_rownames('X1') %>%
    as.matrix()

  save(targets, counts, file = selected_data_file)
}else { load(selected_data_file) }


janitor::tabyl(targets, Group) %>% kable()
Group n percent
ELI 10 0.3333333
NIT 10 0.3333333
SFI 10 0.3333333

Tabla 1: Muestras por grupo

targets %>% datatable()

Tabla 2: Muestras seleccionadas

Llegados a este punto vemos que tenemos cargados el número correcto de muestras por grupo

Filtrado inicial

Antes de comenzar a explorar los datos vamos a eliminar los conteos muy bajos. Eliminaremos todos los genes que a lo largo de todas las muestras tengan menos de 10 observaciones. Aunque el sistema es un poco simple, la idea es simplemente reducir la carga computacional y eliminar el ruido más elemental, a la hora de analizar los datos, DSeq se encargará de eliminar los conteos más ruidosos y además utilizaremos log fold shrinkage para reportar los resultados y penalizar los datos más ruidosos.

filter_lc <- function(mat, threshold = 1) {
  mat[rowSums(mat) > threshold,]
}

counts <- filter_lc(counts, 10)

Control de calidad y exploración

Vamos en primer lugar a controlar la calidad de los datos y a realizar una pequeña exploración de los mismos.

Antes de empezar vamos a cambiar la matriz de conteos y vamos a transformarla en formato largo, deberíamos de ser precavidos ya que esto funcionará bien en un dataset pequeño pero en dataset más grandes podría dar problemas.

counts_long <- counts %>%
  as_tibble() %>%
  rownames_to_column(var = 'gene_id') %>%
  pivot_longer(-gene_id, names_to = 'Sample_Name', values_to = 'count') %>%
  inner_join(dplyr::select(targets, Sample_Name, Group), by = 'Sample_Name') %>%
  mutate(Group = factor(Group, levels = c('NIT', 'SFI', 'ELI'))) %>%
  arrange(Group)

Distribución de conteos por muestras y grupos.

Vamos a comprobar la distribución de conteos por muestra y grupo, primero estudiando su densidad y luego en forma de boxplot.

cowplot::plot_grid(
  ggplot(counts_long, aes(x = log10(count + 1), color = Sample_Name, group = Sample_Name)) +
    geom_density(alpha = .2) +
    lmft +
    theme(legend.position = 'None') +
    labs(title = 'Conteo por muestra'),
  ggplot(counts_long, aes(x = log10(count + 1), color = Sample_Name, group = Sample_Name)) +
    geom_density(alpha = .2) +
    facet_wrap(. ~ Group) +
    lmft +
    theme(legend.position = 'None') +
    labs(title = 'Conteo por grupos'),
  nrow = 2
)

Figura 1: Plots de densidad de los valores de conteo por muestra y grupo

ggplot(counts_long, aes(y = log2(count + 1), x = glue::glue('{Group}_{Sample_Name}'), fill = Group, Group = Group)) +
  geom_boxplot(alpha = .2) +
  lmft +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(title = 'Conteo por muestra', x = 'Muestra')

Figura 2: Boxplot de los valores de conteo por muestra y grupo

Vamos ahora a controlar el conteo total por cada muestra

ggplot(counts_long %>%
         mutate(Sample_Name = glue::glue('{Group}_{Sample_Name}')) %>%
         group_by(Sample_Name) %>%
         summarise(count = sum(count)),
       aes(x = Sample_Name, y = (count / 1e6))) +
  geom_bar(stat = 'identity', alpha = 0, color = 'black') +
  lmft +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

Figura 3: Histograma del número total de conteos por por muestra y grupo

Parece que tenemos una muestra del grupo ELI que está ligeramente desplazada respecto de las otras muestras y además tiene muchas menos medidas. En concreto es la muestra que termina en 2TC6J.

Como no tenemos acceso a los datos iniciales no podemos comprobar cual ha sido el fallo en la muestra, por tanto lo que vamos a hacer es eliminarla y obtener una muestra nueva de nuestra base de datos.

# Remove wrong sample
wrong_sample <- targets[str_detect(targets[['Sample_Name']], '2TC6J'), 'Sample_Name']
targets <- targets[!str_detect(targets[['Sample_Name']], '2TC6J'),]

# Replace sample
new_sample <- read_csv('../data/targets.csv') %>%
  dplyr::filter(Group == 'ELI', Sample_Name != wrong_sample, !(Sample_Name %in% targets[['Sample_Name']])) %>%
  sample_n(1)

targets <- bind_rows(targets, new_sample) %>%
  mutate(rn = Sample_Name) %>%
  column_to_rownames('rn')

# Recreate counts
counts <- read_delim('../data/counts.csv', delim = ';') %>%
  dplyr::select(X1, rownames(targets)) %>%
  column_to_rownames('X1') %>%
  as.matrix()

# Refilter
counts <- filter_lc(counts, 1)

Vamos a volver a hacer el control de calidad con histogramas

counts_long <- counts %>%
  as_tibble() %>%
  rownames_to_column(var = 'gene_id') %>%
  pivot_longer(-gene_id, names_to = 'Sample_Name', values_to = 'count') %>%
  inner_join(dplyr::select(targets, Sample_Name, Group), by = 'Sample_Name') %>%
  mutate(Group = factor(Group, levels = c('NIT', 'SFI', 'ELI'))) %>%
  arrange(Group)
cowplot::plot_grid(
  ggplot(counts_long, aes(x = log10(count + 1), color = Sample_Name, group = Sample_Name)) +
    geom_density(alpha = .2) +
    lmft +
    theme(legend.position = 'None') +
    labs(title = 'Conteo por muestra'),
  ggplot(counts_long, aes(x = log10(count + 1), color = Sample_Name, group = Sample_Name)) +
    geom_density(alpha = .2) +
    facet_wrap(. ~ Group) +
    lmft +
    theme(legend.position = 'None') +
    labs(title = 'Conteo por grupos'),
  nrow = 1
)

Figura 4: Plot de densidad de valor de conteos por muestra y grupo (Muestra Corregida)

ggplot(counts_long, aes(y = log2(count + 1), x = glue::glue('{Group}_{Sample_Name}'), fill = Group, Group = Group)) +
  geom_boxplot(alpha = .2) +
  lmft +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(title = 'Conteo por muestra', x = 'Muestra')

Figura 5: Boxplot de conteos por muestra y grupo (Muestra Corregida)

ggplot(counts_long %>%
         mutate(Sample_Name = glue::glue('{Group}_{Sample_Name}')) %>%
         group_by(Sample_Name) %>%
         summarise(count = sum(count)),
       aes(x = Sample_Name, y = (count / 1e6))) +
  geom_bar(stat = 'identity', alpha = 0, color = 'black') +
  lmft +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

Figura 5: Boxplot de conteos por muestra y grupo (Muestra Corregida)

Parece que hemos corregido el problema que teníamos con nuestra muestra, por tanto proseguimos con el análisis. De ahora en adelante siempre utilizaremos la muestra corregida.

DSeqDataset

Llegado a este punto vamos a transformar el set de datos en un DSeqDataset. Antes de realizar este paso vamos a transformar el campo Group, el de interés para nuestro análisis, en un factor cuyo primer nivel será NIT que será el de referencia cuando analicemos los datos. Lo haremos en este punto para poder transformar los datos en el siguiente paso.

En este punto vamos a añadir los GRanges al objeto DSeq, es importante recordar que la base de datos correcta para realizar esta anotación es la v75, que es la que contiene los genes que aparecen analizados aquí.

targets[['Group']] <- factor(targets[['Group']], levels = c('NIT', 'SFI', 'ELI'))
dds <- DESeqDataSetFromMatrix(countData = counts, colData = targets, design = ~Group)
names_genes_dds <- str_replace_all(names(rowRanges(dds)), '\\..*', '')
genes_ensemb <- genes(EnsDb.Hsapiens.v75)
rowRanges(dds) <- genes_ensemb[names_genes_dds,]

Transformación

Algunos de los métodos de exploración que vamos a utilizar funcionan mejor bajo condiciones de homocedasticidad. Sin embargo esto no es completamente cierto para nuestros datos ya que la varianza aumenta en función del conteo por gen, es decir los genes con mayor conteo tendrán mayor variabilidad.

Podemos verlo claramente en nuestros datos.

ggplot(counts_long %>%
         group_by(gene_id) %>%
         summarise(sd_gene = sd(count), rank_gene = mean(count)) %>%
         mutate(rank_gene = row_number(rank_gene))
  , aes(x = rank_gene, y = sd_gene)) +
  geom_point() +
  ylim(0, 2^10)

Figura 6: Scatter plot entre los genes ordenados de acuerdo a su conteo y su desviación standard. Nótese que el gráfico está cortado en su parte superior.

Como vemos en la Figura 6la variabilidad aumenta con la media.

Vamos por tanto a transformar nuestros datos para reducir la heterocedasticidad. Como son pocas muestras utilizaremos variance stabilizing transformation o vst. Lo haremos con la opción blind que evitará que en el computo de la transformación influya la varianza debida a las condiciones experimentales.

vsd <- vst(dds, blind = F)

Vamos a comprobar que hemos arreglado este problema y la variabilidad de los datos no depende, o lo hace menos de la media.

ggplot(assay(vsd) %>%
         as_tibble(rownames = 'gene_id') %>%
         pivot_longer(-gene_id, values_to = 'count', names_to = 'Sample_Name') %>%
         group_by(gene_id) %>%
         summarise(sd_gene = sd(count), rank_gene = mean(count)) %>%
         mutate(rank_gene = row_number(rank_gene))
  , aes(x = rank_gene, y = sd_gene)) +
  geom_point()

Figura 7: Scatter plot entre los genes ordenados de acuerdo a su conteo y su desviación standard. Nótese que el gráfico al contrario que el anterior, no está cortado en su parte superior.

Como vemos en la Figura 7, aunque no es completamente perfecto, la variabilidad esá menos asociada al conteo por gen. Vamos ahora a explorar los datos.

Distancia entre muestras

sampleDists <- dist(t(assay(vsd)))
library("pheatmap")
library("RColorBrewer")

sampleDistMatrix <- as.matrix(sampleDists)
rownames(sampleDistMatrix) <- vsd$Group
colnames(sampleDistMatrix) <- NULL
colors <- colorRampPalette(rev(brewer.pal(9, "Blues")))(255)
pheatmap(sampleDistMatrix,
         clustering_distance_rows = sampleDists,
         clustering_distance_cols = sampleDists,
         col = colors)

Figura 8: Plot de distancia entre genes un clustering jerárquico.

El agrupamiento es más fuerte entre las muestras del grupo ELI mientras que el resto no parecen agruparse particularmente bien.

PCA

plotPCA(vsd, intgroup = 'Group')

Figura 9: Análisis de componentes principales. Los colores representan los diferentes grupos experimentales.

En el caso del PCA parece más prometedor, existe una diferencia clara entre condiciones en el eje X y además la cantidad de varianza explicada por ese factor es elevada, un 59%.

MDS

mds <- as.data.frame(colData(vsd)) %>%
  cbind(cmdscale(sampleDistMatrix))
ggplot(mds, aes(x = `1`, y = `2`, color = Group)) +
  geom_point(size = 3) +
  coord_fixed()

Figura 10: Plot de Multidimensional Scaling.

Como en el caso de las distancias parece que el grupo ELI se diferencia claramente de los otros mientras que SFI y NIT se mezclan más. Esto es coherente ya que Multidimensional scaling intenta situar los puntos en el plano de tal manera que tengan una distancia igual a la matriz de distancias que hemos calculado previamente.

Agrupamiento de genes

Vamos a seleccionar los 20 genes más variables a través de las muestras de los datos transformados sobre vst

topVarGenes <- head(order(rowVars(assay(vsd)), decreasing = TRUE), 20)

mat <- assay(vsd)[topVarGenes,]
mat <- mat - rowMeans(mat)
anno <- as.data.frame(colData(vsd)[, c("Group", 'Group')])
pheatmap(mat, annotation_col = anno)

Figura 11: Plot de variabilidad entre los 20 genes con clustering jerárquico.

El único cluster que parece asociado a las condiciones está en la esquina superior izquierda. Parece que contiene un set de genes asociados a la condición ELI. EL cluster central no parece particularmente asociado a ninguna condición.

Conclusiones

Parece que los datos son correctos y la exploración de los datos también indica que podemos proceder con el análisis de expresión diferencial.

Análisis de expresión diferencial (DSeq) y Anotación

Como vamos a aplicar DSeq no es necesario aplicar ningún tipo de transformación por lo tanto trabajaremos directemente sobre los conteos, en el objeto DSeqDataset que hemos creado en el apartado anterior. En este caso no hemos añadido la variable sex al modelo, sin embargo sería correcto comprobar contra un modelo con y sin esta variable para saber si es necesario introducirla en el modelo, o quizás incluirla de todos modos ya que la influencia del sexo es probablemente cierta, se podría consultar con el investigador. Como la PEC no indica

Si quisiésemos por ejemplo trabajar con limma deberíamos de transformar los datos (VOOM) ya que no está preparada para conteos.

En este caso disponemos de tres condiciones:

  • NTI: Non-Infiltrated Tissues
  • SFI: Small Focal Infiltrates
  • ELI: Extensive Lymphoid Infiltrates

Vamos a aplicar varios enfoques en este caso:

Compararemos las dos condiciones experimentales (SFI y ELI) contra la condición control (NTI) y luego compararemos las condiciones experimentales entre ellas.

En todo los casos utilizaremos como corte de significación p<.05 y en los casos de FDR un 5%. En otro caso se indicará explícitamente.

dds_pair_file <- 'dds_pair.RData'
if (!file.exists(dds_pair_file)) {
  dds_pair <- DESeq(dds, parallel = T)
  save(dds_pair, file = dds_pair_file) }else {
  load(dds_pair_file)
}

Extraemos los resultados para cada contraste de interés y luego los plotearemos uno por uno. También vamos a calcular los resultados con un shrinkage de el fold que nos ayudará a corregir la variabilidad en los casos en el que conteo sea muy bajo o muy variable y podremos ordenar mejor los resultados.

En el caso de los contrastes contra el nivel de referencia utilizaremos apeglm para el shrinkage y ashr, para el caso del ELI y SFI.

Para realizar las anotaciones de los genes utilizaremos una vez mas la Ensembl v75 en lugar de la org.Hs.eg.db que hemos utilizado en la PEC anterior.

annot <- function(rs) {
  if(class(rs)=='GRanges'){
    k <- rownames(mcols(rs))
  }else{
    k <- rownames(rs)
  }

  rs$symbol <- mapIds(EnsDb.Hsapiens.v75,
                      keys = str_replace_all(k, '\\..*', ''),
                      column = "SYMBOL",
                      keytype = "GENEID",
                      multiVals = "first")
  rs$entrez <- mapIds(EnsDb.Hsapiens.v75,
                      keys = str_replace_all(k, '\\..*', ''),
                      column = "ENTREZID",
                      keytype = "GENEID",
                      multiVals = "first")
  return(rs)
}

rs_dds_pair_file <- 'rs_dds_pair.RData'
if (!file.exists(rs_dds_pair_file)) {


  k <-rs_SFI_NIT <- annot(results(dds_pair, name = 'Group_SFI_vs_NIT', alpha = .05))
  rs_ELI_NIT <- annot(results(dds_pair, name = 'Group_ELI_vs_NIT', alpha = .05))
  rs_SFI_ELI <- annot(results(dds_pair, contrast = c('Group', 'SFI', 'ELI'), alpha = .05))

  rs_SFI_NIT_LFC <- annot(lfcShrink(dds_pair, coef = 'Group_SFI_vs_NIT', type = 'apeglm'))
  rs_ELI_NIT_LFC <- annot(lfcShrink(dds_pair, coef = 'Group_ELI_vs_NIT', type = 'apeglm'))
  rs_SFI_ELI_LFC <- annot(lfcShrink(dds_pair, contrast = c('Group', 'SFI', 'ELI'), alpha = .05, type = 'ashr'))

  rs <- bind_rows(mutate(as_tibble(rs_SFI_NIT, rownames = 'gene'), contrast = 'Group_SFI_vs_NIT'),
                  mutate(as_tibble(rs_ELI_NIT, rownames = 'gene'), contrast = 'Group_ELI_vs_NIT'),
                  mutate(as_tibble(rs_SFI_ELI, rownames = 'gene'), contrast = 'Group_SFI_vs_ELI'))

  rs_LFC <- bind_rows(mutate(as_tibble(rs_SFI_NIT_LFC, rownames = 'gene'), contrast = 'Group_SFI_vs_NIT'),
                      mutate(as_tibble(rs_ELI_NIT_LFC, rownames = 'gene'), contrast = 'Group_ELI_vs_NIT'),
                      mutate(as_tibble(rs_SFI_ELI_LFC, rownames = 'gene'), contrast = 'Group_SFI_vs_ELI'))
  save(rs_SFI_NIT, rs_ELI_NIT, rs_SFI_ELI,
       rs_SFI_NIT_LFC, rs_ELI_NIT_LFC, rs_SFI_ELI_LFC,
       rs, rs_LFC,
       file = rs_dds_pair_file) }else {
  load(rs_dds_pair_file)
}

Hemos calculado también los resultados con log fold shrinkage para observar como determinados genes con grandes cambios en el log fold van asociados a bajos conteos y alta variabilidad. LFCShrink penaliza estos genes reduciendo su log fold. Utilizaremos estos últimos para reportar los resultados. Además también hemos anotado todos los genes

Resultados

Annotations

Como hemos visto que hay muchos genes que no tienen anotación en org.Hs.eg, vamos primero a hacer un conteo rápido para saber cuantos genes mapeados tienen un correspondiente Symbol o EntrezID.

rs_SFI_NIT %>%
  as_tibble(rownames = 'gene_id') %>%
  transmute('Symbol' = ifelse(!is.na(rs_SFI_NIT[['symbol']])|rs_SFI_NIT[['symbol']]=='', 'Yes', 'No'),
            'Symbol' = factor(Symbol, levels = c('Yes', 'No')),
            'Entrez' = ifelse(!is.na(rs_SFI_NIT[['entrez']]), 'Yes', 'No'),
            'Entrez' = factor(Entrez, levels = c('Yes', 'No'))
  ) %>%
  janitor::tabyl(Symbol, Entrez) %>%
  janitor::adorn_totals(where = c('row', 'col')) %>%
  janitor::adorn_title()
        Entrez            
 Symbol    Yes    No Total
    Yes  22762 20649 43411
     No      0     0     0
  Total  22762 20649 43411

Tabla 3: Tabla de éxito de mapeado de Ensembl a EntrezID y Symbol

Como vemos aproximadamente la mitad de los genes no tienen ‘EntrezID’ de acuerdo a la base de datos mientras que todos tienen un Symbol asociado.

Tablas de resultados

Vamos a presentar los resultados en tablas que son ordenables y filtrables para cada uno de los contrastes por separado. Hemos eliminado las filas cuyos p.valores correspondían a NA. Los datos de cada contraste se dividen en Log Fold Positivo y Negativo.

table_report <- function(df) {
  one_table <- function(df, fun) {
    df %>%
      as_tibble(rownames = 'gene_id') %>%
      fun() %>%
      arrange(padj) %>%
      drop_na(-symbol, -entrez) %>%
      sample_n(100) %>%
      mutate_if(is.numeric, function(x) { signif(x, digits = 3) })
  }

  return(list(one_table(df, fun = function(x) { dplyr::filter(x, log2FoldChange > 0) }),
              one_table(df, fun = function(x) { dplyr::filter(x, log2FoldChange < 0) })))

}

short_sum <- function(df, alpha = .05) {
  df %>%
    as_tibble(rownames = 'gene_id') %>%
    transmute(Regulated = if_else(log2FoldChange > 0, 'Up', 'Down'),
              Regulated = factor(Regulated, levels = c('Up', 'Down')),
              Significant = if_else(!is.na(padj),
                                    if_else(padj < alpha, 'Yes', 'No'),
                                    'Low Count'),
              Significant = factor(Significant, levels = c('Yes', 'No', 'Low Count'))
    ) %>%
    janitor::tabyl(Regulated, Significant) %>%
    janitor::adorn_totals(where = c('row', 'col')) %>%
    janitor::adorn_title()
}

print_datatable <- function(x, caption) { datatable(x %>% arrange(padj), filter = 'top', caption = caption) }

Tabla SFI vs. NIT

short_sum(rs_SFI_NIT_LFC, .05) %>% kable()
Significant
Regulated Yes No Low Count Total
Up 181 15330 7403 22914
Down 3 13585 6909 20497
Total 184 28915 14312 43411 
t <- table_report(rs_SFI_NIT_LFC)
print_datatable(t[[1]], caption = 'Log Fold Positivo')
print_datatable(t[[2]], caption = 'Log Fold Negativo')

Tabla ELI vs. NIT

short_sum(rs_ELI_NIT_LFC, .05) %>% kable()
Significant
Regulated Yes No Low Count Total
Up 3403 13882 6643 23928
Down 1258 13081 5144 19483
Total 4661 26963 11787 43411 
t <- table_report(rs_ELI_NIT_LFC)
print_datatable(t[[1]], caption = 'Log Fold Positivo')
print_datatable(t[[2]], caption = 'Log Fold Negativo')

Tabla SFI vs. ELI

short_sum(rs_SFI_ELI_LFC) %>% kable()
Significant
Regulated Yes No Low Count Total
Up 1207 12892 6093 20192
Down 2817 13025 7377 23219
Total 4024 25917 13470 43411 
t <- table_report(rs_SFI_ELI_LFC)
print_datatable(t[[1]], caption = 'Log Fold Positivo')
print_datatable(t[[2]], caption = 'Log Fold Negativo')

Plots

Volcano plot

my_volcano <- function(df, alpha = .05, lim_fold = 1, title = '') {

  df <- df %>%
    tidylog::drop_na() %>%
    mutate(sig_p = pvalue < alpha, sig_adj_p = padj < alpha) %>%
    tidylog::drop_na()

  unadj_p <- ggplot(df, aes(x = log2FoldChange, y = -log10(pvalue), color = sig_p)) +
    geom_point() +
    geom_hline(yintercept = -log10(alpha), linetype = 2) +
    geom_vline(xintercept = lim_fold) +
    geom_vline(xintercept = -lim_fold) +
    labs(x = 'LogFold',
         y = '-log10(p)',
         title = paste(title, 'Sin corregir por FDR'),
         color = paste0('p<', alpha)) +
    facet_wrap(. ~ contrast) +
    theme(legend.position = 'none') +
    lmft

  adj_p <- ggplot(df, aes(x = log2FoldChange, y = -log10(padj), color = sig_adj_p)) +
    geom_point() +
    ggthemes::theme_tufte(base_family = 'sans') +
    geom_hline(yintercept = -log10(alpha), linetype = 2) +
    geom_vline(xintercept = lim_fold) +
    geom_vline(xintercept = -lim_fold) +
    labs(x = 'LogFold',
         y = '-log10(p_adj)',
         color = paste0('p<', alpha),
         title = paste(title, 'Corregido por FDR p-value')
    ) +
    facet_wrap(. ~ contrast) +
    lmft

  return(list(unadj_plot = unadj_p, adj_plot = adj_p))
}
cowplot::plot_grid(plotlist = my_volcano(df = rs, alpha = .05, lim_fold = 2, title = 'Sin LFC'), nrow = 2)

cowplot::plot_grid(plotlist = my_volcano(df = rs_LFC, alpha = .05, lim_fold = 2, title = 'Con LFC'), nrow = 2)

fig_nums(name = 'plot_SFI_NIT', caption = 'Volcano SFI vs. NTI con y sin LFC', display = F)

Figura 13: Volcano plot para cada una de las condiciones y los resultados con y sin shrinkage. Las líneas verticakes indican un logfold de 2 y la horizontal el límite de significación al 5%

La línea horizontal indica el límite de significación para un 5%, las líneas verticales un fold de 2, los colores indican si un gen ha resultado significativo o no.

Como vemos en el caso sin shrinkage aparecen multitud de genes con un fold muy elevado, sin embargo al aplicar LFC el número se reduce drásticamente, probablemente sea indicativo de abundantes genes con variabilidad elevada. La única comparación que no sale fuertemente penalizada es SFI vs. ELI. En la sección anterior hemos visto los volcano plots que hemos utilizado en primer lugar para comparar la diferencia entre los resultados cuando utilizamos un shrinkage y cuando no lo utilizamos. Por tanto no los repetiremos aquí

MA-Plot

Hemos ploteado para todos los contrastes tanto los resultados sin como con shrinkage. En el caso que no hemos un shrinkage podemos comprobar como la estimación del log fold es más variable con conteos bajos (logFold altos pero sin embargo no son significativos) y sin embargo se vuelve menos variable cuando aumenta el conteo. Una consecuencia de esto es que logFold de una determinada magnitud no son significativos si el conteo es bajo y sin embargo si lo son con conteos alto. En rojo aparecen coloreados los resultados significativos. Vemos también como la transformación vst es efectiva ya que reduce el logFold en los conteos bajos, con lo que en esa franja de los gráficos se concentran en cero.

SFI vs. NIT

plotMA(rs_SFI_NIT, ylim = c(-5, 5))

plotMA(rs_SFI_NIT_LFC, ylim = c(-5, 5))

Figura 14: MA plot de los resultados con y sin shrinkage

ELI vs. NIT

plotMA(rs_ELI_NIT, ylim = c(-5, 5))

plotMA(rs_ELI_NIT_LFC, ylim = c(-5, 5))

Figura 15: MA plot de los resultados con y sin shrinkage

SFI vs. ELI

plotMA(rs_SFI_ELI, ylim = c(-5, 5))

plotMA(rs_SFI_ELI_LFC, ylim = c(-5, 5))

Figura 15: MA plot de los resultados con y sin shrinkage

Plot en el espacio genómico

Como ejemplo de gráfico en el espacio genómico vamos a mostrar el gen más significativo en la comparación SFI vs. ELI. Para ello primero reanalizamos los datos pero lo hacemos en formato GRanges y sacamos los fold changes de los resultados con shrinkage.

Vamos a plotear las 100,000 base alrededor del gen más significativo.

Utilizaremos tracks para mostrar: - Los p.valores ajustados (-log10) - Los log2fold corregidos - Los genes significativos al 5% - Los genes cuyo log2fold absoluto corregido es mayor que 2

resGR <- annot(results(dds_pair, name = "Group_SFI_vs_NIT", format = "GRanges"))
resGR$log2FoldChange <- rs_SFI_ELI_LFC$log2FoldChange
resGR$logpadj <- -log10(resGR$padj)

window <- trim(resGR[which.min(resGR$padj)] + 1e5)
strand(window) <- '*'
resGRsub <- resGR[resGR %over% window]
naOrDup <- is.na(resGRsub$symbol) | duplicated(resGRsub$symbol)
resGRsub$group <- ifelse(naOrDup, names(resGRsub), resGRsub$symbol)

status_p <- factor(ifelse(resGRsub$padj < 0.05 & !is.na(resGRsub$padj), "sig", "notsig"))
status_2 <- factor(ifelse(abs(resGRsub$log2FoldChange) > 2, "sig", "notsig"))

options(ucscChromosomeNames = FALSE)
g <- GenomeAxisTrack()
p_annot <- AnnotationTrack(resGRsub, name = "significant - gene ranges (padj<.05)", feature = status_p)
log_annot <- AnnotationTrack(resGRsub, name = "logfold > 2 - gene ranges", feature = status_2)
p_data <- DataTrack(resGRsub, data = "logpadj", baseline = 0,
               type = "h", name = "-log10(padj)", strand = "+")
log_data <- DataTrack(resGRsub, data = "log2FoldChange", baseline = 0,
               type = "p", name = "log2 fold change", strand = "+")

plotTracks(list(g, p_data, log_data, p_annot, log_annot), groupAnnotation = "group",
           notsig = "deepskyblue", sig = "tan2")

Figura 16: Plot track del gen más significativo y las 100,000 bases vecinas. En narajan aparecen los genes signficativos al 5% (track superior) o con un log fold superio a 2 (track inferior).

Podemos observar como en este área la mayoría de los genes tienen un log fold inferior a cero por tanto su conteo es menor en la condición ELI que en la condición SFI y existen varios genes significativos alrededor. También podemos observar como incluso genes con un fold change elevado no son significativos, probablemente debido a la variabilidad.

Análisis de enriquecimiento

  • Hemos utilizado el paquete clusterProfiler

    • Para estos análisis hemos seleccionado todos aquellos genes cuyo p-valor corregido con FDR es menor a .1
    • Hemos realizado un análisis por contraste
    • Hemos utilizado las bases de datos KEGG y GO
    • En el caso de KEGG hemos corregido siempre con un p<.05 y un q<.05 para ORA y p<.05 para GSEA
    • En el caso de GO hemos corregido siempre con un p<.05 y un q<.05, y p<.05 para GSEA.
    • Hemos realizado 10000 permutaciones en los GSE

    Enlaces a los tipos de gráficos utilizados:

En este sección utilizaremos dos tipos de análisis sobre dos bases de datos diferentes:

  • Overrepresentation y Gene Set Enrichment Analysis
    • Gene Ontology GO y Kyoto Encyclopedia of Genes and Genomes, KEGG

Interpretación de los resultados

Para los resultados de cada uno de estos análisis presentamos el mismo conjunto de tablas y resultados. Una tabla que muestra los resultados obtenidos, con información relevante sobre cada uno de ellos y un enlace a cada uno de los elementos de la GO o KEGG encontrados. Los tipos de gráficos varían según el tipo de análisis.

En el caso de KEGG cuando una ruta resulta significativa, el diagrama de la ruta se descarga y se colorea dependiendo de la expresión de los genes incluidos en la ruta para el análisis.

enrich_set <- rs_LFC %>%
  drop_na() %>%
  dplyr::select(p.fdr = padj, ENTREZID = entrez, fold = log2FoldChange, contrast) %>%
        dplyr::filter(p.fdr < .05)

to_entrezid <- function(df) {
  df %>% dplyr::pull(ENTREZID)
}

to_fold <- function(df) {
  sort(tibble::deframe(df %>% dplyr::select(ENTREZID, fold)), decreasing = T)
}

enrich_kegg <- function(cont, gene_list, pvaluecutoff = .05, qvaluecutoff = .05) {
  gene_list_contrast <- to_entrezid(gene_list %>% dplyr::filter(contrast == cont))
  gene_list_kegg <- bitr_kegg(gene_list_contrast,
                              fromType = 'ncbi-geneid',
                              toType = 'kegg', organism = 'hsa')$kegg
  enrich_res <- enrichKEGG(gene_list_kegg,
                           organism = 'hsa',
                           pvalueCutoff = pvaluecutoff,
                           keyType = 'kegg',
                           qvalueCutoff = qvaluecutoff)
  # Ok con el ncbi-geneid! https://www.genome.jp/dbget-bin/www_bget?dme:Dmel_CG14941 Ejemplo
  return(enrich_res)
}

cont <- 'Group_SFI_vs_NIT'
gene_list <- enrich_set

gse_kegg <- function(cont, gene_list, pvaluecutoff = .05, n_perm = 10000) {
  gene_list_contrast <- to_fold(gene_list %>% dplyr::filter(contrast == cont))
  names(gene_list_contrast) <- bitr_kegg(names(gene_list_contrast),
                                         fromType = 'ncbi-geneid',
                                         toType = 'kegg', organism = 'hsa')$kegg
  enrich_res <- tryCatch({ gseKEGG(gene_list_contrast,
                                   organism = 'hsa',
                                   pvalueCutoff = pvaluecutoff,
                                   nPerm = n_perm,
                                   keyType = 'kegg') },
                         error = function(cond) { tibble(ID = character(),
                                                         Description = character(),
                                                         SetSize = numeric(),
                                                         enrichmentScore = numeric(),
                                                         NES = numeric(),
                                                         pvalue = numeric(),
                                                         p.adjsut = numeric(),
                                                         qvalues = numeric(),
                                                         rank = numeric(),
                                                         leading_edge = character()
                         ) })
  # Ok con el ncbi-geneid! https://www.genome.jp/dbget-bin/www_bget?dme:Dmel_CG14941 Ejemplo
  return(enrich_res)
}

enrich_go <- function(cont, gene_list, pvaluecutoff = .05, qvaluecutoff = .05) {
  gene_list_contrast <- to_entrezid(gene_list %>% dplyr::filter(contrast == cont))
  enrich_res <- enrichGO(gene = gene_list_contrast,
                         OrgDb = org.Hs.eg.db,
                         keyType = 'ENTREZID',
                         ont = "CC",
                         pAdjustMethod = "BH",
                         pvalueCutoff = pvaluecutoff,
                         qvalueCutoff = qvaluecutoff)
  enrich_res <- setReadable(enrich_res, OrgDb = org.Hs.eg.db, keyType = "ENTREZID")
  return(enrich_res)
}

cont <- 'Group_SFI_vs_NIT'

gse_go <- function(cont, gene_list, n_perm = 10000, pvaluecutoff = 0.05) {
  gene_list_contrast <- to_fold(gene_list %>% dplyr::filter(contrast == cont))
  enrich_res <- tryCatch({ gseGO(geneList = gene_list_contrast,
                                 OrgDb = org.Hs.eg.db,
                                 keyType = 'ENTREZID',
                                 ont = "CC",
                                 nPerm = n_perm,
                                 minGSSize = 100,
                                 maxGSSize = 500,
                                 pvalueCutoff = pvaluecutoff,
                                 verbose = FALSE) },
                         error = function(cond) { tibble(ID = character(),
                                                         Description = character(),
                                                         GeneRatio = character(),
                                                         BgRatio = character(),
                                                         pvalue = numeric(),
                                                         p.adjust = numeric(),
                                                         q.value = numeric(),
                                                         Count = numeric()) }
  )
  if (nrow(enrich_res) > 0) {
    enrich_res <- setReadable(enrich_res, OrgDb = org.Hs.eg.db, keyType = "ENTREZID")
  }
  return(enrich_res)
}
plot_enrich_ORA <- function(enrich_df, contrast, type) {
  if (nrow(enrich_df) == 0) {
    return(invisible(NA))
  }
  p1 <- heatplot(enrich_df) + labs(title = paste(contrast))
  p2 <- dotplot(enrich_df)

  plot(p1)
  plot(p2)

  if (type == 'GO' && nrow(enrich_go_SFI) > 1) {
    p3 <- goplot(enrich_df)
    plot(p3)
    # cowplot::plot_grid(p1,p2,p3, labels = 'AUTO', nrow = 3)
  }
}

plot_enrich_GSE <- function(enrich_df, contrast, type) {
  if (nrow(enrich_df) == 0) {
    return(invisible(NA))
  }

  if (type == 'GO') {
    p1 <- emapplot(enrich_df) + labs(title = paste(contrast))
  }else {
    p1 <- ridgeplot(enrich_df) + labs(title = paste(contrast))
  }
  p2 <- dotplot(enrich_df) + labs(title = paste(contrast))
  # cowplot::plot_grid(p1, p2, labels = 'AUTO', nrow=2)
  plot(p1)
  plot(p2)
}

enrich_go_table <- function(enrich_go) {
  suppressWarnings(enrich_go %>%
                     tibble::as_tibble() %>%
                     dplyr::mutate(ID = glue::glue("<a href='http://amigo.geneontology.org/amigo/term/{ID}' target='_blank'>{ID}</a>", ID = ID)) %>%
                     dplyr::mutate_if(is.numeric, signif, digits = 2) %>%
                     dplyr::select(-one_of("geneID", "core_enrichment")) %>%
                     DT::datatable(escape = FALSE))
}


enrich_kegg_table <- function(enrich_kegg) {

  if (nrow(enrich_kegg) != 0) {

    enrich_kegg_tibble <- enrich_kegg %>% tibble::as_tibble()

    gen_url <- function(x, pathID) { paste0("http://www.kegg.jp/kegg-bin/show_pathway?",
                                            pathID, "/", x[pathID, "geneID"]) }

    suppressWarnings(enrich_kegg_tibble %>%
                       dplyr::mutate(Link = lapply(ID, function(x) { stringr::str_sub(gen_url(x = enrich_kegg, pathID = x), 1, -2) }),
                                     Description = paste0("<a href='", Link, "' target='_blank'>", Description, "</a>")) %>%
                       dplyr::mutate_if(is.numeric, signif, digits = 2) %>%
                       dplyr::select(-Link, -one_of("geneID", "core_enrichment")) %>%
                       DT::datatable(escape = FALSE)) }else
    { tibble() %>% DT::datatable(escape = FALSE) }
}

enrich_kegg_pathview <- function(enrich_kegg, gene_list, dir_path = './kegg_pathview') {

  enrich_kegg_tibble <- enrich_kegg %>% tibble::as_tibble()

  if (nrow(enrich_kegg_tibble) == 0) { return(invisible(NULL)) }

  dev.null <- lapply(enrich_kegg_tibble$ID, function(x) { if (!file.exists(file.path(dir_path, paste0(x, '.xml')))) {
    pathview(gene.data = to_fold(gene_list),
             pathway.id = x,
             species = 'hsa',
             kegg.dir = dir_path
    )
    file.copy(from = paste0(x, ".pathview.png"),
              to = file.path(dir_path, paste0(x, ".pathview.png")))
    file.remove(paste0(x, ".pathview.png"))

  } })

  dev.null <- lapply(enrich_kegg_tibble$ID, function(x) { plot(cowplot::ggdraw() + cowplot::draw_image(image_trim(image_read(file.path(dir_path, paste0(x, ".pathview.png")))))) })

}

Gene ontology

Over-representation analysis

enrich_go_SFI <- suppressMessages(suppressWarnings(enrich_go('Group_SFI_vs_NIT', enrich_set)))
enrich_go_ELI <- suppressMessages(suppressWarnings(enrich_go('Group_ELI_vs_NIT', enrich_set)))
enrich_go_SFI_ELI <- suppressMessages(suppressWarnings(enrich_go('Group_SFI_vs_ELI', enrich_set)))

SFI vs. NIT

enrich_go_table(enrich_go_SFI)
plot_enrich_ORA(enrich_go_SFI, 'Group_SFI_vs_NIT', type = 'GO')

#### ELI vs. NIT

enrich_go_table(enrich_go_ELI)
plot_enrich_ORA(enrich_go_ELI, 'Group_ELI_vs_NIT', type = 'GO')

SFI vs. ELI

enrich_go_table(enrich_go_SFI_ELI)
plot_enrich_ORA(enrich_go_SFI_ELI, 'Group_SFI_vs_ELI', type = 'GO')

Gene Set enrichment analysis

gse_go_SFI_ELI <- suppressMessages(suppressWarnings(gse_go('Group_SIF_vs_ELI', enrich_set)))
gse_go_ELI <- suppressMessages(suppressWarnings(gse_go('Group_ELI_vs_NIT', enrich_set)))
gse_go_SFI <- suppressMessages(suppressWarnings(gse_go('Group_SFI_vs_NIT', enrich_set)))

SFI vs. NIT

enrich_go_table(gse_go_SFI)
plot_enrich_GSE(gse_go_SFI, 'Group_SFI_vs_NIT', type = 'GO')

ELI vs. NIT

enrich_go_table(gse_go_ELI)
plot_enrich_GSE(gse_go_ELI, 'Group_ELI_vs_NIT', type = 'GO')
Resultados de GSE para GO

Resultados de GSE para GO

Resultados de GSE para GO

Resultados de GSE para GO

SFI vs. ELI

enrich_go_table(gse_go_SFI_ELI)
plot_enrich_GSE(gse_go_SFI_ELI, 'Group_SFI_vs_ELI', type = 'GO')

KEGG

En el caso de KEGG se han ocultado los pathways yo que en algunos casos son demasiados y hacen dificil navegar el documento. Se pueden desplegar en cada sección haciendo click en los puntos indicados.

Over-representation analysis

enrich_kegg_SFI <- suppressMessages(suppressWarnings(enrich_kegg('Group_SFI_vs_NIT', enrich_set)))
enrich_kegg_ELI <- suppressMessages(suppressWarnings(enrich_kegg('Group_ELI_vs_NIT', enrich_set)))
enrich_kegg_SFI_ELI <- suppressMessages(suppressWarnings(enrich_kegg('Group_SFI_vs_ELI', enrich_set)))

SFI vs. NIT

enrich_kegg_table(enrich_kegg_SFI)

Haga click para ver los pathways 

enrich_kegg_pathview(enrich_kegg_SFI, enrich_set)
plot_enrich_ORA(enrich_kegg_SFI, 'Group_SFI_vs_NIT', type = 'KEGG')

ELI vs. NIT

enrich_kegg_table(enrich_kegg_ELI)

Haga click para ver los pathways 

enrich_kegg_pathview(enrich_kegg_ELI, enrich_set)
plot_enrich_ORA(enrich_kegg_ELI, 'Group_ELI_vs_NIT', type = 'KEGG')

SFI vs. ELI

enrich_kegg_table(enrich_kegg_SFI_ELI)

Haga click para ver los pathways 

enrich_kegg_pathview(enrich_kegg_SFI_ELI, enrich_set)

plot_enrich_ORA(enrich_kegg_SFI_ELI, 'Group_SFI_vs_ELI', type = 'KEGG')

Gene Set enrichment analysis

gse_kegg_SFI_ELI <- suppressMessages(suppressWarnings(gse_kegg('Group_SIF_vs_ELI', enrich_set)))
gse_kegg_ELI <- suppressMessages(suppressWarnings(gse_kegg('Group_ELI_vs_NIT', enrich_set)))
gse_kegg_SFI <- suppressMessages(suppressWarnings(gse_kegg('Group_SFI_vs_NIT', enrich_set)))

SFI vs. NIT

enrich_kegg_table(gse_kegg_SFI)

Haga click para ver los pathways 

enrich_kegg_pathview(gse_kegg_SFI, enrich_set)
plot_enrich_GSE(gse_kegg_SFI, 'Group_SFI_vs_NIT', type = 'KEGG')

ELI vs. NIT

enrich_kegg_table(gse_kegg_ELI)

Haga click para ver los pathways 

enrich_kegg_pathview(gse_kegg_ELI, enrich_set)
plot_enrich_GSE(gse_kegg_ELI, 'Group_ELI_vs_NIT', type = 'KEGG')

SFI vs. ELI

enrich_kegg_table(gse_kegg_SFI_ELI)

Haga click para ver los pathways 

enrich_kegg_pathview(gse_kegg_SFI_ELI, enrich_set)
plot_enrich_GSE(gse_kegg_SFI_ELI, 'Group_SFI_vs_ELI', type = 'KEGG')

Efectos ocultos

En este caso tenemos una variable oculta bastante evidente que es el sexo, vamos a comprobar si somos capaces de encontrarla con SVA

library("sva")

dat  <- counts(dds_pair, normalized = T)
idx  <- rowMeans(dat) > 1
dat  <- dat[idx, ]
mod  <- model.matrix(~ Group, colData(dds_pair))
mod0 <- model.matrix(~   1, colData(dds_pair))
svseq <- svaseq(dat, mod, mod0, n.sv = 2)
Number of significant surrogate variables is:  2 
Iteration (out of 5 ):1  2  3  4  5  
par(mfrow = c(2, 1), mar = c(3,5,3,1))
for (i in 1:2) {
  stripchart(svseq$sv[, i] ~ dds_pair$sex, vertical = TRUE, main = paste0("SV", i))
  abline(h = 0)
}

Figura 17: Plot de las variables ocultas, en el eje x la variable real sex que no estaba incluida en el modelo.

No parece que la separación sea completamente perfecta, pero si parece que hay una separación aproximada particularmente en la primera variable oculta.

Archivos generados

Los archivos aparecen comprimidos en el repositorio en formato zip.

save_dose_result = function(enrich_df, fname){
  enrich_df %>% tibble::as_tibble() %>% readr::write_csv(x = ., path = paste0(fname,'.csv'))
}

if(!file.exists('../results.zip')){

  # Tables
  rs_SFI_NIT %>% as_tibble(rownames = 'gene_id') %>% write_csv('../results/results_SFI_NIT.csv')
  rs_ELI_NIT %>% as_tibble(rownames = 'gene_id') %>% write_csv('../results/results_ELI_NIT.csv')
  rs_SFI_ELI %>% as_tibble(rownames = 'gene_id') %>% write_csv('../results/results_SFI_ELI.csv')

  rs_SFI_NIT_LFC %>% as_tibble(rownames = 'gene_id') %>% write_csv('../results/results_SFI_NIT_LFC.csv')
  rs_ELI_NIT_LFC %>% as_tibble(rownames = 'gene_id') %>% write_csv('../results/results_ELI_NIT_LFC.csv')
  rs_SFI_ELI_LFC %>% as_tibble(rownames = 'gene_id') %>% write_csv('../results/results_SFI_ELI_LFC.csv')

  # Enrichement
  enrich_go_SFI %>% save_dose_result(fname='../results/ora_go_SFI')
  enrich_go_ELI %>% save_dose_result(fname='../results/ora_go_ELI')
  enrich_go_SFI_ELI %>% save_dose_result(fname='../results/ora_go_SFI_ELI')

  gse_go_SFI %>% save_dose_result(fname='../results/gse_go_SFI_ELI_ELI')
  gse_go_ELI %>% save_dose_result(fname='../results/gse_go_ELI')
  gse_go_SFI_ELI %>% save_dose_result(fname='../results/gse_go_SFI_ELI')

  enrich_kegg_SFI  %>% save_dose_result(fname='../results/ora_kegg_SFI')
  enrich_kegg_ELI %>% save_dose_result(fname='../results/ora_kegg_ELI')
  enrich_kegg_SFI_ELI %>% save_dose_result(fname='../results/ora_kegg_SFI_ELI')

  gse_kegg_SFI  %>% save_dose_result(fname='../results/gse_kegg_SFI')
  gse_kegg_ELI %>% save_dose_result(fname='../results/gse_kegg_ELI')
  gse_kegg_SFI_ELI %>% save_dose_result(fname='../results/gse_kegg_SFI_ELI')

  zip::zipr(files='../results',zipfile = '../results.zip')
}
Nombre
results_SFI_NIT.csv
results_ELI_NIT.csv
results_SFI_ELI.csv
results_SFI_NIT_LFC.csv
results_ELI_NIT_LFC.csv
results_SFI_ELI_LFC.csv
ora_go_SFI.csv
ora_go_ELI.csv
ora_go_SFI_ELI.csv
gse_go_SFI_ELI_ELI.csv
gse_go_ELI.csv
gse_go_SFI_ELI.csv
ora_kegg_SFI.csv
ora_kegg_ELI.csv
ora_kegg_SFI_ELI.csv
gse_kegg_SFI.csv
gse_kegg_ELI.csv
gse_kegg_SFI_ELI.csv

Apéndice de código

Todo el código generado aparece en el markdown y se puede consultar o ocultar